Why Do Granular Materials Stiffen with Shear Rate? A Test of Novel Stress-Based Statistics 
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Recent experiments exhibit a rate-dependence for granular shear such that the stress grows linearly in the 
logarithm of the shear rate, y. Assuming a generalized activated process mechanism, we show that these obser- 
vations are consistent with a recent proposal for a stress-based statistical ensemble. By contrast, predictions for 
rate-dependence using conventional energy-based statistical mechanics to describe activated processes, predicts 
a rate dependence that of (ln(y)) 1 / 2 . 



Understanding disordered solids, such as foams, glasses, 
polymers, colloids and granular materials is a great challenge 
for statistical physics. Several of these systems, including 
granular materials, fall outside the rubric of conventional sta- 
tistical mechanics because they are dissipative. But, these ma- 
terials exhibit well defined statistical distributions. Several 
novel approaches have been recently proposed(TJ O [3] |U [5j 
13 [SI to characterize the statistics of these dissipative ma- 
terials,. We focus on testing one of the proposed statistical 
frameworks, the force or stress-based ensemblesO 171 0, 
and specifically the stressed-based ensemble hypothesized by 
SH and BC(7] [8) to account for the coupling between forces 
and geometry. Here, we test this hypothesis by showing that it 
can account for experimentally observed logarithmic strength- 
ening with increasing shear rate in slowly sheared granular 
systems. Many models of slow, dense granular flows assume 
that the internal stresses are independent of shear rate. Lin- 
ear rate-dependence for shear stresses occurs for Newtonian 
fluids. Non-linear dependence on rate is common in non- 
Newtonian fluids, as well as in the "glassy" systems noted 
aboveO [TO). The stress-based ensemble offers an explicit 
framework for analyzing the rheology of such non-thermal, 
glassy systems. 

Recent 2D [11] and 3D lfT2l experiments on sheared dense 
granular materials, showed mean stresses that grow linearly 
with ln(y), where y is the shear rate. Indeed, rate-dependence 
spans many decades in y, as seen in Fig.[T] which show time- 
averaged stresses acquired in a 2D Couette shear experiment 
that is sketched in Fig. [I] Fig. [2] shows typical traces of stress 
vs. time for several different shear rates, y. These data have 
also been acquired for various packing fractions, (]), relatively 
near the critical packing fraction, ty C9 below which the system 
is unj ammed fT3l [T7 l and shearing stops. 

We construct a model for the behavior of the force fluctua- 
tions seen in the Couette shear experiments. The first model 
premise is that the shearing process causes build-up of inho- 
mogeneous stress structures, such as force chains, which fail 
when they reach a critical yield stress. A blow-up of a typical 
time series for the stress in Fig. [3] gives a good flavor of this 
process. 

The actual failure process is complex except very near jam- 
ming (<|> ~ <|) c ), where there are typically only one or two vis- 



ible force chains in an observation region. More generally, 
for § much larger than there is a (strong) force network, 
and when a stress drop occurs, it is typically concentrated on 
a subset of grains on a more or less linear region of the force 
network, often extending for many grains. The stress drop is 
very fast relative to the build-up, indicating a failure of part 
of the network, while much of the remaining network is only 
weakly affected. We use the common term "force chains" to 
refer to segments of the strong force network that exhibit the 
force build-up and failure process that is the fundamental ori- 
gin of granular force fluctuations, but we do not consider the 
specific origins of the failures (e.g. shear-transformation-zone 
events fl"5lL force chain buckling fT6lL etc.). Rather, we explore 
the role played by the background of force fluctuations gener- 
ated during the build-up and failure process. We ask whether 
the effect of these fluctuations can be described in terms of 
a stress-based ensemble. And we exploit the fact that force 
changes following failure are localized to force chains, with a 
much weaker effect on the rest of the system. 

We expect that a chain will fail if the force/stress on it, 
a, exceeds a characteristic value G m . This is reminiscent of 
Coulomb failure, but refers to failures of localized structures, 
with no strict frictional analogy. We begin by considering a 
single event consisting of the birth-to-death cycle of a single 
force chain, relevant to systems with § ~§ c , and then return 
to an accounting for multiple events occuring in a given ob- 
servation region for packings with § » 4> c . 

A key premise of the model is that the failure is an acti- 
vated process aided by stress fluctuations, similar to a ther- 
mally activated escape from a potential well. The potential 
well is replaced by a stress trap that models the meso-scale 
strong force network, and chain failures correspond to escape 
from a trap. The fluctuations of thermal equilibrium are re- 
placed by fluctuations of stress in the network, characterized 
by a temperature-like quantity |'8'|. These fluctuations occur as 
the granular assembly moves through a series of states at me- 
chanical equilibrium. For thermally activated processes, the 
rate of escape is proportional to e~$ Ebarrier , where P = (ksT)~ l 
is the inverse temperature and iWner is the barrier height. 
To construct a framework for activated dynamics in systems 
where the fluctuations are athermal stress fluctuations, we ap- 
peal to a recently-developed statistical framework for granular 
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FIG. 1: Top, sketch of Couette shear apparatus, not to scale. 
R = 19.2cm. Data is for small/medium shearing wheels, where 
Ri = 6.7 cm/5.3 cm. Particles are either circular or pentagonal in 
cross section, and have dimensions between 0.7 cm and 0.9 cm. Bot- 
tom: Data for the mean stress in a segment of a 2D granular Couette 
experiment containing roughly 200 particles vs. shear y. Data (nor- 
malized by Gdc — 4.11Nm _1 ) are for different packing fractions (]) 
relative to the critical packing fraction, ty c , where stesses fall to zero. 
Symbols are: circles: § — § c = 0.0012, squares: § — § c = 0.0091, dia- 
monds: — § c = 0.0152, upward pointing triangles: cf) — 4> c = 0.0189, 
left-pointing triangles: circles: § — § c = 0.0226. We also consider 
other data (below) from experiments by Hartley and Behringer| 11 1. 
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FIG. 2: Time series for different y (in mHz) for the stress, for the 
smaller Ri above. 



FIG. 3: Blow up of time series data for the stress at y = 0.0533mHz 
showing stress fluctuations over 0.2 rotations of the inner shearing 
wheel. The data indicate the approximately linear loading up of force 
structures and their more rapid failure. 



assemblies(7][8]]. In this framework, the ensemble of mechan- 
ically stable states is defined by a Boltzman-like probability 
distribution: 



P v = (l/Z)e- aT ^ 



(1) 



where T v = S v c v is an extensive quantity related to the stress 
of the configuration v, and S v is the area occupied by the 
grains QH). In Eq. [T] a is the analog of the inverse temper- 
ature, (3, and characterizes the fluctuations. In analogy with 
thermally activated processes, the probability per unit time of 
chain failure is then given by 



Pf=Aexp[-ar barrier ]. 



(2) 



Assuming that the area fluctuations are small compared to the 
stress fluctuations (the area is fixed in the experiments dis- 
cussed above), the effective barrier to be surmounted by a 
force chain with a stress a on it becomes Tb arr i er = S(o m — a), 
and the failure rate (per unit time) is: 



Pf=Aexp[-(G m -a)/o )]. 



(3) 



We have used Go to denote 1 / (aS). 

Both A, the attempt frequency, and G may depend on y, but 
to lowest order, we will treat these as constants. We expect 
that the stress on a force chain increases linearly in time until 
a chain fails (e.g. Fig. [5]). Thus, 



c = a't = Hyt, 



(4) 



where £ is a measure of rate of stress increase per unit shear 
deformation, which we also assume is a constant. With this 
picture, the force chain loads up steadily in time, but the 
probability of failure depends on the closeness of a to o m . 
If a << a m , the probability of failure/unit- time should be 
low, but as a approaches o m , the probability of failure should 
become large. A range of parameters where the process is 
strongly activated is A « y, and a m >> o Q , the analog of 
the low-temperature limit of a thermally activated process. 
This limit likely applies to the experiments since the time- 
dependence of the stress is dominated by slow build up and 
rapid release. This is the regime we focus on in this work. 
The assumption of strongly activated behavior is also consis- 
tent with the fact that the experiments suggest logarithmic rate 
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dependence over many decades of y, without a crossover to 
different behavior at low y 

We treat the loading up of the networks in a probabilistic 
manner. If the probability for the chain to survive until time t 
without failing is P s (t), then 

P s {t + dt)=P s (t){l-P f (t)dt), (5) 

or 

P~ l dP s /dt = -P f . (6) 

We first consider the idealized situation of isolated chains, so 
that we need only consider one such chain at any given instant, 
and focus on the evolution of stresses from the formation of 
the chain to its failure. The time average moments of G are: 

d* = (i/r) [ T dt(a(t))« 

Jo 

= J d%P s (x)P f (x)J\t(c(t)) n / 1 dzP s (z)P f (z)z(l) 

where we have used, T = f dxP s (x)Pf(x)x. Defining (a n ) = 
fdzP s (z)P f (z)(a(x)) n , we can write G n = (o n + l )/((n + 
1 ) (a) ) The above equation reflects an ensemble average over 
chains surviving up to time x and failing within x to x + dx, 
with probability P s (x)Pf(x). We can use Eq.[6]to write: 

(a n ) = -J\a(x)r(dP s /dx)dx (8) 

Introducing the dimensionless time, = yu, we can write 

d\n(P s )/dQ = -£exp(Le/G ), (9) 

where we we have absorbed some of the constants into a sin- 
gle expression B = (A/y) exp(— c m /o ). Integrating Eq. [9] 
using P s (0) = 1, and integrating Eq.[8]by parts, we can write 
(o n ) as: 

(o n )=nZ n J d n - l exp[(B/p)(l-exp(pd))]dd (10) 

Here we define p = £/g . In the special case n = 1, it is 
possible to relate this integral to known functions ifTTl . but this 
does not appear to be true for the general case, and it is now 
convenient to introduce a dimensionless rate, s = p/B: 

s = (yE)exp (c m /c )/(Ac ) (11) 

Differentiating (g h ) with respect to s yields, after a bit of al- 
gebra: 

sd { G n ) jds = - {a n ) /s + na { G n ~ 1 ) . (1 2) 

To calculate the time-averaged stress a, we need to calculate 
(a) and (a 2 ). The average (a), 

d(G)/ds=-(c)/s 2 + c /s. (13) 



We are concerned with the strongly activated regime of A < < 
y and a m >> Go, s » 1, which justifies an expansion involv- 
ing s and its logs. In lowest order when s « 1: 

<G)^G ln(s)+C (14) 

or 

(g) =G m + G ln(E/G )+G ln(y/A)+Ci, (15) 

where C\ is a constant of integration. Similarly, keeping the 
leading terms in s, the expression for (g 2 ) is: 

d(a 2 )/ds ~ 2g (o)/s 

(g 2 ) = G 2 (\ns) 2 + C 2 (16) 

where C2 is another integration constant. Assuming the inte- 
gration constants are much smaller than hx(s)) 9 

(g 2 ) 

G=^-G \ns (17) 

In fact, we can now see if omitting the term — (o) / s 2 at lowest 
order is self consistent. We estimate (g) ~G m , and then obtain 
a ratio of the two terms on the right side of Eq. [13] 

[(a)/^]/\L/(sp)] ~ (G m /Z)(A/y)exp(-G m /G ). (18) 

Indeed, this ratio should be small, so the dominant rate effect 
on the mean stress should be a logarithmic stengthening. 

We now turn to how the mean-i.e. long-time averaged 
stresses from the whole network within a measurement do- 
main would be manifested in a continuous shear experi- 
ment. For instance, in the 2D Couette experiments [11], time- 
dependent force data are obtained in a finite region compris- 
ing roughly 10% of the whole system. In 3D experiments lfT2l . 
pressure measurements are made over an area which contacts 
several tens of particles. In both cases, the mean stress is com- 
puted as a time average. Here, we focus on the 2D experiment, 
since it has yielded data over a range of densities. If (|) were 
such that on average only one force chain existed at a time, 
then our analysis so far would yield the mean stress. How- 
ever, in general, it is necessary to adjust this result for the 
mean number of chains (or failure events), N(§), that are gen- 
erated per unit time interval or angular displacement in the 
region of interest. Note that to lowest order, this is simply 
a function of (|). In principle, we can determine this quan- 
tity and hence G in order to compare to preditions from the 
stress ensemble For instance, from recent experiments 

by Sperl et al. on the same 2D Couette svstem fT8lL we es- 
timate (at fixed y) that N(ty) ((|> — § c ) a and G °c (((|) — § c ) b , 
where a ^ 2 and b ~ 1. Data from Hartley et al. fffl . such as 
Fig.[T] yield N(§)g by taking the slope of g/g^c vs. ln(y). 
Unfortunately, the current data is not sufficiently precise to 
give a good determination of G by this process. 

Nevetheless, Eq.18, with the above ^-dependent correc- 
tions, makes a prediction, the key result of this work: G oc 
ln(y), with a proportionality constant of G , multiplied by 
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FIG. 4: Fits (solid lines) to Eq. [l9]of experimental data for Ac at (|> — 
§c = 0.01373, and y's of 0.066 mHz (circle), 0.6645 mHz (square) 
and 13.307 mHz (triangle). The fits yield Go ~ 2.0, independent of y 
and o° /o ~ 1.58 ,1.82 ,2.0 (lowest to the highest value of y.) 



N(§). Qualitatively, the strengthening of the material with 
shear rate occurs because the failure probability per unit time 
for a given (a m — g)/g q is independent of y whereas the 
buildup of stress grows linearly with y. There is, therefore, 
a competition between the stress buildup and the failure pro- 
cess. If the failure rate were independent of a (large Go 
limit), then, it follows from Eq. ^ that the average stress 
would grow linearly with the shear rate, as in Newtonian flu- 
ids. In the opposite limit of failure occurring only if the force 
chains are loaded up to G m (Go — > limit), there would be 
a fixed time °c 1/y between failures and the average stress 
would become independent of the shear rate. The exponen- 
tial increase of the failure rate with g leads to the logarith- 
mic strengthening. Interestingly, Eq. [12] indicates a much 
weaker rate dependence for the variance, V = (g 2 ) — ((c)) 2 : 
dV/ds = (-V + ((a)) 2 )/s 2 . 



The distribution of the stress drops, Aa, occuring for each 
"avalanche" in a time series is sensitive to the distribution of 
G m , and this distribution can be calculated exactly ||T0j [T9) for 
an experimentally realistic (20) exponential distribution of G m : 
<?-o m /o? n . For Zy/Go << A(e Aa / a ° - 1), the distribution is 



P(Ao) 



1 



i (l+ao/O 



e Ao/o _ i - 



(19) 



which we have fitted to the experimental data for stress drops, 
ensuring that the chosen data meet the criterion for which Eq. 
19 is applicable. Representative data at — § c = 0.01373 
are shown in Fig. [4] The ratio of the two fitting parame- 
ters, Go/g^, follows by fitting the exponent of the power law 
regime, and G^ can be estimated from a fit to the exponen- 
tial tail. We use these values to constrain the fitting to the full 
form of Eq. [19] 

The analysis presented in this letter is reminiscent of the 
soft glassy rheology (SGR) model(9] [TO) with the important 
distinction that the role of energy is being played by the stress. 
The noise temperature x of the SGR is replaced by 1 /a in the 
stress-based ensemble. The SGR model incorporates disorder 
through a distribution of activation barriers, which in the con- 
text of the current work would translate to a distribution of G m . 



The average stress is not sensitive to the distribution of G m in 
the large s limit, and scales as (my). The difference with the 
SGR result, (xlny) 1 / 2 , can be explained through replacement 
of an energy barrier by a stress barrier! 19 ]. 

The present analysis is a first step towards understanding 
linear logarithmic strengthening in granular materials which 
distinguishes granular rheology, with dissipative interactions, 
from that of other materials. A key point is that a stress ensem- 
ble rather than a energy ensemble yields the correct rate scal- 
ing. Force networks are visually obvious in the experiments, 
but their quantitative connection to the complete stress states 
has not yet established. Characterization of mesoscale struc- 
tures such as force networks, and their connection to macro- 
scopic variables such as stress remains a great challenge for 
the field. Adopting the framework of the SGR model with its 
meanfield approach to correlations but with the noise tempera- 
ture and energy replaced by their counterparts from the stress- 
based ensemble should provide a fruitful avenue for building 
these connections. 
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